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The three-site Bose-Hubbard model subject to atom losses: the boson-pair dissipation 

channel and failure of the mean-field approach 
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We employ the perturbation series expansion for derivation of the reduced master equations for 
the three-site Bose-Hubbard model subject to strong atom losses from the central site. The model 
describes a condensate trapped in a triple-well potential subject to externally controlled removal of 
atoms. We find that the 7r-phase state of the coherent superposition between the side wells decays via 
two dissipation channels, the single-boson channel (similar to the externally applied dissipation) and 
the boson-pair channel. The quantum derivation is compared to the classical adiabatic elimination 
within the mean-field approximation. We find that the boson-pair dissipation channel is not captured 
by the mean-field model, whereas the single-boson channel is described by it. Moreover, there is a 
matching condition between the zero-point energy bias of the side wells and the nonlinear interaction 
parameter which separates the regions where either the single-boson or the boson-pair dissipation 
channel dominate. Our results indicate that the M-site Bose-Hubbard models, for M > 2, subject 
to atom losses may require an analysis which goes beyond the usual mean-field approximation for 
correct description of their dissipative features. This is an important result in view of the recent 
experimental works on the single site addressability of condensates trapped in optical lattices. 

PACS numbers: 03.75.Gg; 03.75.Lm; 03.75.Nt 



I. INTRODUCTION 

The mean-field approximation, formally obtained by 
replacing the boson creation and annihilation operators 
by complex scalars, is usually employed for description of 
bosonic many-body systems when the number of bosons 
is large (for instance, see Refs. [ll-Q)- Such a replace- 
ment can be justified by the \/]V-scaling of the boson 
operators, that is, when the populations are large, the 
commutators - the source of quantum corrections - are 
negligible |l|. The relation between the quantum and 
mean-field descriptions is a subject of intensive studies. 
The quantum description is necessary for the bifurca- 
tions, which modify significantly the quantum spectrum 
[J] , the quantum collapses and revivals [5[ , and the many- 
body quantum corrections to the mean- field theory (y, |7j . 
Making explicit the iV-scaling of the operators and iden- 
tifying the iV-scaling of the parameters for a fixed particle 
density, reveals the link of the mean-field approximation 
to the Wentzel-Kramers-Brillouin semiclassical approach 
to the discrete Schrodinger equation J8J, now employed 
in the Fock space with the inverse number of bosons 1/N 
playing the role of an effective Planck constant (see, for 
example, Ref. Q). Therefore, the mean-field limit, as the 
semiclassical limit of a discrete Schrodinger equation, is 
also singular. Hence, besides the pronounced quantum 
corrections/fluctuations at the bifurcations/instabilities, 
one must be prepared to find even a qualitative disagree- 
ment between the mean-field description and the full 
quantum consideration even when the populations are 
large, as it is the case, for instance, in the nonlinear boson 
model which describes tunneling of boson pairs between 
two modes, see Refs. (lOllllj. 

The main purpose of the present paper is to study 



the dissipation dynamics of the atom filled lattice sites 
coupled to a common dissipated sites. Our motivation 
is the recent advancement in the quantum microscopy 
techniques and the current experiments on the single 
site addressability in the optical lattices [Tig, LU1 , where 
controlled atom losses are induced in selected sites of 
a two-dimensional optical lattice. We develop a direct 
method based on the perturbation expansion for deriva- 
tion of the reduced quantum master equations for the 
Bose-Hubbard models with dissipation (we consider the 
atom losses) and compare the quantum derivation with 
the mean-field description. We concentrate on the three- 
site Bose-Hubbard model, which is the simplest model 
describing atom filled sites coupled to a common dis- 
sipated one and describes, for instance, cold atoms or 
Bose-Einstein condensate trapped in a triple-well poten- 
tial subject to removal of atoms from the central well, 
see Fig. [T] The three-site Bose-Hubbard model was also 
noted to possess many features of complexity of a general 
quantum dynamics, as the Wigner-Dyson spectral statis- 
tics and quantum chaos [1J, [l5[ ■ It is also the simplest 
model where the boson-pair tunneling, originating from 
the nonlincarity of the model, is possible. 

We derive the reduced quantum master equations for 
the coherent modes describing the condensate in the side 
wells, then consider the mean-field approach and com- 
pare the results of the two approaches. We note here 
that we consider an open quantum system and, as such, 
described by the quantum master equation [161 Il7j ] . How- 
ever, in the case of the atom losses, the mean- field formu- 
lation is straightforward (this also applies to the case of 
the m ultip le-site atom losses of the experimental works of 
Refs. [lj, O)- Contrary to the fact that the mean-field 
approximation applies with a good accuracy to the two- 
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FIG. 1. (Color online) The three-site Bose-Hubbard model 
setup. A laser/electron beam removes atoms from the central 
well with a rate V. The ej represents the zero-point energy 
of the jih well. The two quantum channels of dissipation of 
the 7r-phase coherent mode, described in the text, are shown 
schematically by arrows. 



site Bose-Hubbard model with atoms losses and a noise 
[la . KL9| , it is shown here that the correct and complete 
description of the three-site model (in general, the M-site 
Bose-Hubbard models, with M > 3) requires going be- 
yond the usual mean-field approach. This disagreement 
stems from the quantum boson-pair dissipation channel, 
due to the nonlinear interaction (this is similar to the 
boson-pair tunneling resulting in the qualitative failure of 
the mean- field approach in Ref. [HI]). Moreover, there is 
a matching condition between the zero-point energy bias 
of the side wells and the nonlinear interaction parame- 
ter which separates the regions where cither the single- 
boson or the boson-pair dissipation channel dominate. 
Hence, one has to use the full quantum consideration, i.e. 
the quantum master equation reduction methods, to de- 
scribe the decay of the subsystem (in our case, the quan- 
tum modes describing the side wells) , which then may be 
treated with further approximations, even resembling the 
mean-field approach. However, the point is that without 
invoking the quantum consideration at some stage, i.e. 
working just within the mean-field approach, one will be 
unable to describe the dissipativc behavior of the filled 
sites coupled to a common dissipated site, which con- 
clusion is also relevant to the recent experiments of Ref. 

The paper is organized as follows. In section |TT] we 
formulate the quantum master equation. The derivation 
of the reduced master equation for the side modes (i.e., 
the coherent superposition modes over the left and right 
wells in Fig. [IJ is given in section IIIII In section IIVI a 
similar reduction is applied to the reduced master equa- 
tion of section IIIII producing the master equation for the 
7r- phase coherent mode. In section IVl the adiabatic elim- 
ination within the mean-field approximation is studied. 
The concluding remarks and discussion is contained in 



II. THE PROBLEM FORMULATION IN TERMS 
OF THE MASTER EQUATION 



A quantum channel representation for a local atom 
losses (e.g., from a single lattice site ), s ee Fig. [T] can 
be given in the Fock space as follows [l9[ 

\kj)\0)B -> Vp\kj - 1)|1)h + V^^P\k 3 )\0)B, 

\kj) is the ket- vector of the Fock space of a dissipat- 
ing boson mode, \X)r describes the atom counter, and 
p = p(kj,t) is the probability. Note that p(kj,5t), for 
small St, depends linearly on the number of atoms in the 
dissipating mode. As the result, a Lindblad term with 
the generator L = yTcij appears in the master equa- 
tion for the density matrix. Here T is the dissipation 
rate parameter. We consider the dissipation acting on 
the central well (denoted with j = 3 in Fig. [IJ of the 
triple- well trap, thus the master equation for the density 
matrix reads 



dp »rrr 1 

Tt = —h [H ' p] 



1 



a 3P% - 2 



n 3 p 



\pn, 



(1) 



By H we denote the three-site Bose-Hubbard Hamilto- 
nian, 



H = -J{[a\+a\]a z +al[a 1 +a 2 \)+^e :j a]a J +Uy^\a]fa 



\»j) a j, 



3=1 j=l 

(2) 
where a, and aj (j = 1,2,3) are the local boson modes 
describing occupation of the respective well, J is the lin- 
ear tunneling rate, Sj is the zero-point energy of the re- 
spective well, and U is the atomic interaction parameter 
proportional to the s-wave scattering length. 

Due to the linear coupling of the central well to the 
side wells in the Bose-Hubbard model, it is convenient 
to use the new canonical basis b\ = (ai + a 2 )/\/2, 
b-2 = (ai — a 2 )/\/2, and 63 = a 3 . Here the modes bi }2 
describe, respectively, the zero-phase and 7r-phase coher- 
ent superpositions between the side wells. The Bose- 
Hubbard Hamiltonian becomes 

■Jis(b% 

^^([m + n 2 \ 2 + [b\b 2 



H 



blh) ~ J 12 (b% 



b\bi) 

bth} 2 



en 3 



(3) 



where n,- 



6'-6i and we have introduced the parameters 
(e 2 - ei)/2, e = e 3 - (£1 + £ 2 )/2 



J13 = v2J, J12 

(and dropped an inessential term proportional to N 

ni +n 2 +TI3). 



III. THE REDUCED MASTER EQUATION FOR 
THE SIDE WELLS 



In the strongly dissipative case, T >• J13/S, the central 
well is emptied on the time scale t ~ 1/T, while the co- 



herent modes 6i_and b 2 almost retain their populations 
(see also Ref. [19||). It turns out that in this case one 
can derive an approximate master equation for the side 
wells alone valid for the times t S> 1/r, when the central- 
well dissipation has the fastest time-scale in the system. 
Introducing the small parameter e = J\ 3 /hT we require 
that J 12 /hT = 0(e) and U(m t2 )/hT = 0(e). This ap- 
proximation is valid for an arbitrary bias e between the 
central and side wells. We will see that, already in the 
first order in e, the density matrix for the whole system 
is not factorized and the central- well population strongly 
depends on that of the side wells. Nevertheless, in the 
higher orders of the perturbation theory (namely, in the 
second and third orders in e), the central well acts as an 
effective reservoir for the side wells leading to the master 
equation in the standard Lindblad form. 

First of all, let us pass to the interaction picture. In- 
troduce 

Hj{t) = U$(t)HiU (t), p = U^(t)p(t)U (t) 

U Q (t) = U 12 (t) ® U 3 (t) = exp{- l -{H 12 +H 3 )At}, (4) 

where At = t — to, to is some initial time, and H\ 2l H 3 
and Hi are parts of the Hamiltonian ^ describing the 
side wells, the central well and the interaction between 
them, respectively (we have also subtracted the nonessen- 
tial term U(ni+n 2 +n 3 ) from the system Hamiltonian to 
simplify the presentation). The master equation in the 
interaction picture reads 



renormalization correction to the first term. The coef- 
ficient c\ has the following meaning: Tr 3 {6 3 (t)p 3 (t)} = 
dfi(t)e~~ At + 0(e~ rAt ), with /i being the eigenvalue 
of the unitary transformation in Eq. §5§: C/3 (t) 1 1 ) = 
fi(t)|l), in our case /i(t) = e ~ ieAt/h (note also that 

63 (t) = /i(*)6a). 

Inserting Eq. © into Eq. |J7J| we get the general solu- 
tion in the following form (suggested by the form of p 3 
itself) 



P W (t)=V(t)\p^ 



F t (*) + 0(e-5 At ), (10) 



.(0+1) _ Jo) , m 



where p 12 — p\ 2 + p\ 2 and V(t) reads 

V(t) = fflcp{oi(t)6i(t)6j(t)} = 7+a 1 (t)& 1 (t)6 3 (t)+0(e 2 ), 

(11) 
with a scalar function ct\{t) = 0(e). Indeed, taking 
derivative of the solution (fTTj) and using Eqs. (|7|) and 
©, we obtain the equation to be satisfied 

^ = ^>|0>(0|yt + ^ 2 )^|0)(0|^ 



dt dt 



(1) 



dt 



+v^®mw"= ~[ fl >(*)^«« 



TD{Vpf 2 +1) ®\Q) (O|W} + 0(e 



(12) 



$ = -Urn), P }+T (k(t)pbi(t) - Wp \ 



dt 



:P n 3 



where b 3 (t) = [/g(t)& 3 [/o(t). Below we will work in the 
interaction picture and drop the hat, for simplicity. The 
density matrix expansion reads p 



Pi° 2 } « Pf 



P 



(i) 



pi 2 ) _|_ 0(e 3 ), where, for instance, p\ 2 = Tr 3 {p} (an arbi- 
trary non-factorized expansion, p^ = J^.. c ijP\ 2 i ® Paj * 
leads to the same result). In the lowest orders in e we 
have: 



dPri _ dp { ° ] 



dt 



dt 



TD(t){p^}, 



(6) 



We evaluate 



(o) 



-j:[Hi(t),p^ 



= ! f i (/i(*)&i(*)Pi° 2 ) ®|l)(0|-/i.c) (13) 

and, using Eq. ([TT]) . 

r^(t){^° 2 +1) ®|0)(0|^ } 

= -| (aiWA'WftiC*)^? ® Il>(0| + ft-c) + 0(e 2 ). (14) 
Eqs. and (TT3l give immediately 



^ - -^[ff/W,^ ® /4 0) ] + rz?(t){p< 1 )}, (7) 



^ 22) =-^Tr 3 {[^(t),pW]}. 



dt 






(0) • 



The solution to the equation for p 3 is as follows 



(8) 



p { 3 \t) = \0}(0\ + e-i At [c 1 \l)(0\+h.c.}+O(e- rAt ), (9) 



dp 12 
dt 



0(e— At ). 



(15) 



Since J 12 /HT = 0(e) and U{m, 2 )/hT = 0(e), the deriva- 
tive of b\{t) is of order e, 



rffri(t) 

dt 



■~{H l2 ,h(t)]=0(e). 



(16) 



Now, by taking into account Eqs. (|TT1) . (fT3|) - (|T6|) . one 
sees that Eq. (fi"2j) is satisfied by setting 



where the last term represents all higher-order exter- 
nal products of the Fock basis states and a decaying 



da 1 
~dt 



r e 

2 +i n 



1 *^ 13 m\ 

Ql + — . (17) 



ion (|17[) gives 




2iJ 13 
ai= hT 


' 2ie 
1 + ^ 



+ 0(e-7 At ). 



(18) 



to those in Eq. ((13)) . Expression ([22)1 also defines the 
general form of p^ 2 \t): 



p^ = B 00 (t) 



B u (i) ®|1)<1| 



The final step is to insert Eq. (fTU|) into Eq. ©, use 
Eqs. ([6]) and (|T5|) and take the trace over the central- well 
subspace, keeping only the terms up to 0{e 2 ). Returning 
to the Schrodinger picture, we get (up to a correction of 
the order C(e 3 )) 



dpi 



dt 



h 

+G(e 



[Hi2,Pl2 
"At 



-[£Rb\h,p\2 



(0)i 



). 



with T R given by 



T R = \ ai \ 2 T = 



h 2 v 



'fir 



(19) 

i 
, (20) 



and Di{p} = bipb\ — ^{b\bi, p}. The interaction- 

induced Lamb shift e R of the zero-point energy of the 
coherent zero-phase mode b\ reads 



e R = -Ji 3 Re(ai) = ~ £ -pr- 



(21) 



Our aim is to further reduce Eq. (|19p to a master equa- 
tion for the coherent 7r-phase mode 62, which has unusual 
dissipation features (see below). To this end, however, 
one has to consider the contribution to Eq. ([T9|) coming 
from the next order in e, i.e. 0(e 3 ). Indeed, similar to 
the derivation of this section, the reduced master equa- 
tion for mode 62 is obtained under the conditions that 
J12, U(ni t 2-) "C hT R , what makes the Hamiltonian part in 
the master equation (|19[) smaller than the Lindblad part, 
hence the former should be discarded in the present order 
C(e 2 ). Thus, in the second-order approximation, the co- 
herent mode 62 has no dissipation dynamics at all (only 
the Hamiltonian evolution described by H2 = ^n 2 ,)- It 
only appears in the higher-order version of the master 
equation (TT9")) , 

To derive the third-order correction to Eq. (|19p . we 
need to find the form of p^ 2 ^(i), which satisfies equation 
similar to Eq. ([7]) , but now with the inhomogeneous term 

-^[^(i),^i 2 +1) ®|0)(0|yt] 



U 



13 



V2a 1 (t)/ 2 *(t)6 2 (t)^° 2 ) ® |2)(0| - h.c. 



2ilm{a 1 (t)}b 1 (t)p ( $ >b\(t) ® |1)(1| 
a^bl^hit)^ ~ h.c.) ®\0)(Q\ 



fi(t)h(t)p$®\l)(0\-h.c. 



0(e 2 



-2»i 



'-At 



(22) 



The 



where U 3 (t)\2) = f 2 (t)\2), i.e. f 2 = e— < — 

first three lines of Eq. (|2"2")l give the terms additional 



+ (B 10 (t) <8> |1)(0| + B 20 (i) <8> |2)(0| + /i.e.) , (23) 

where the operators By(t) = 0(e 2 ) act on the subspace 
of the side wells. They satisfy, in view of Eqs. (|2U1) . (|13l) . 
(fT5| . ([21")) and (|22|l. the following equations: 



rfBn 

dB 00 
dt 



-TBn+TRhi^p^blit), 
rB 11 -~[e R n 1 (t),p ( $] 



,(°) 



,(°). 



■y(ni(t)^+ri>i(<) 

(i) 



Bi = ai(t)/ 1 *(t)6i(tM2. 
£20 = a a (t)/ 2 *(t)&?(t)p£\ 

where c<i(t) is given by Eq. (|18[) and 



d«2 



r e + t/ 
2 + *^T 



1 -V^Jis , . 

a 2 +« — r — ai(i). 



(24a) 

(24b) 

(24c) 
(24d) 

(25) 



One can easily solve Eq. (|24a[) and by using integration 
by parts represent the result as 



r fl, 



B n (t) = -fh(t)p^bl(t) + 0(e 3 ) + 0(e- rAt ). (26) 



Substituting Eq. ([26)) into Eq. (|24b|) and using the mas- 
ter equation (fT^j) one obtains 



(2) 



B 00 (t)=p\y+O(e-^ t ) 
Finally, the solution to Eq. (|25|) reads 



(27) 



2y/2J 2 

+ o( e -5 At ) 



,2e 

'/if 



1 + t 



2(e + f/) 



n -1 



/T 



^l + 0( e 3 ) + ( e -^).( 28 ) 



The fact that the operators B\\ and B20 and B02 do 
not contribute to the equation for p\2 and the explicit 
expressions for B 00 and B 10 , Eqs. ([2"6")l and ([2"T)l , lead to 
the same Lindblad form of the reduced master equation, 
i.e. in the Schrodinger picture we get 

dp 



= - jrpia + £ R ni,pi2] + TrD^pu} + 0(e~ At ), 

(29) 
which is now valid up to correction of order C(e 4 ). Eqs. 
(|23)) . ([26)) -([28 )) also show that the density matrix of the 
full system p can be rewritten as (now in the Schrodinger 
picture) 



p(t) = V[p 12 (t) 



7 T + 0(e d ) + 0(e" = At ), (30) 



where p\ 2 {t) is taken up to the second order in e, V(t) — 
cx.p{aibib\} = I + ot\b\b\ + \ct 2 b\(b\) 2 with a.\ given by 
Eq. flT5J|. Eqs. (J2~S| and (J3DJ) are the main results of this 
section. Obviously, the full density matrix (|50| is not 
factorized, nevertheless, the reduced density matrix of a 
subsystem satisfies the Markovian master equation (|29|) 
in the Lindblad form. Note that the difference between 
the density matrix p\ 2 and the one obtained by tracing 
the full density matrix of Eq. (|30|) is a constant term of 
the second order given by B\\ in Eq. (|26[) . which does 
not contribute to Eq. (f2T))) . Hence, Eq. (|30p is consistent 
with the approximation made. 

In Fig^ [3J we use the Monte Carlo wave-function 

method [20] and find that an excellent agreement of the 
reduced master equation ([2l))) with the full Eq. ([T]) ex- 
tends also to the intermediate values of T (we have there 
e = 0.2). We have also verified that, for large T, the 
modes ai >2 (and, hence, bi t2 ) stay practically unchanged 
while the dissipating mode 63 = 03 is quickly emptied 
(the inset of Fig. [2]). 
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FIG. 2. (Color online) Comparison between the full 3-mode 
Eq. |T| and the reduced 2-mode Eq. (f29|) . The population 
of as is given by the dashed lines. The populations of modes 
a\ and a^ (solid lines) of the 3-mode model are compared to 
that of the 2-mode model (dotted lines). Here F = 5J/ft, 
U — 0.2J, Ej = 0. We have used 5000 quantum trajectories. 
The inset shows the short-time dynamics of populations for 
larger T = 50J/h. The initial state is the Fock state with the 
total of 10 atoms (5 occupying the mode ai and 5 occupying 
the mode 03. 



Let us note the specific features of Eq. ([29]) . We see 
that in the strongly dissipative case, quite similar to Eq. 
([T]), mode b 2 can retain a significant part of its popula- 
tion, while 61 loses almost all atoms, on the time scale 
t ~ l/r^, and after that stays practically empty. In the 
long run, on the time scale much longer that l/r^, the 
population of mode b 2 drops to the single-atom level (this 
is seen also in Fig. [21 where it happens on the short time 
scale due to a small T). 

In the linear case Eq. ([2T?f acts like a dispersive beam- 



splitter (see, for example, Ref. [2l|). Thus, a strong loss 
in the central well induces quantum correlations between 
the side wells, i.e. for times t ^> I/Tr the cold atoms 
occupy the state 



I*, 



$rio>~(ot-<4) n io> 



(31) 



which is unaffected by the dissipation described by Eq. 

In the nonlinear case, the dissipation of the 7r-phasc 
mode 62 is surprisingly non-trivial. This can be clearly 
demonstrated in the case when the decay of the zero- 
phase mode b\ occurs on the faster scale than the inter- 
mode dynamics. To uncover the details, let us use the 
higher-order validity of Eq. (J29|) and reduce it to a 
single-mode master equation for 62, by assuming that 
U(n 2 }/h < T R and J 12 /h < T R . 



THE REDUCED MASTER EQUATION FOR 
tt-PHASE COHERENT MODE 



IV. 



We now reduce the master equation (|29|) to that for 
the mode b 2 alone, in the case of a strong dissipation of 
mode b\ as compared to the Hamiltonian dynamics of 
the coherent modes b\ and b 2 . The derivation is similar 
to that of the previous section. First, we pass to the 
interaction representation: 

H H (t) = uHt)H H U(t), p 12 = U\t) Pl2 {t)U{t) 
U{t) = U x {t) ® U 2 (t) = exp{-^(i/ 1 +iJ 2 )At}, (32) 

where the respective Hamiltonian terms, derived from 
Hamiltonian ([3|) with account of the Lamb shift (|21[) . 
read: 



Hi = ERfl! 



U 



Ho 



U 



H u = -Ji2(b\b 2 + blb 1 ) + 2Un 1 n 2 + ^[b\b 2 ] 2 + ^[blb 1 ] 2 . 

(33) 
Introducing the small parameter e' = 0( J\ 2 /hY r), and 
assuming that U (n 2 ) j (HT #)) = C(e'), one can derive the 
equations for the density matrix in the interaction repre- 
sentation, which turn out to be similar to those of the pre- 
vious section (see Eqs. ©-©) with the obvious changes. 

The form of the density matrix p 12 is also similar to that 
of Eqs. (j7| and pT|) : 

p{\\t) = W(t) [pf +1) ® |0)(0|] W\t) + G(e^ At ), 

(34) 
where W(i) reads 



W(t) 



n(t)b 2 (t)b\(t) + s 2 (t)b 2 2 (t)[b\(t)} 2 , (35) 



with some scalar functions Sj(t) = 0(e'). Using the same 
routine as in the previous section, one obtains the equa- 
tions for the parameters Sj(t): 



ds\ 
~dT ~ 


- 


~Tr 


£r 

h 


, .Jia 


(36a) 


ds 2 

~dt ' ~ ' ~ 


r ^- 2£R 
Fr + i—r- 

h 


U 
S2 - l 2h- 


(36b) 


These can be easily solved to give: 




2zJl2 

nT R 


nT R _ 


-i 


+ 0(e- r -? At ), 





S 2 



iU 



1 



2e R 



R 



2HT R [ HT 

Define the following parameters: 

-l 



0(e 



-i#AtN 



K 2 

3 



h?T 



R 



1 



2e 



R 



hi 



-£r-, 



(37) 



where K\ = 2J\ 2 and K 2 = U. Then, in a similar way 
as in the previous section, by taking the partial trace 
over the subspace of mode b\ , one derives a closed master 
equation for mode 62 . In the Schrodingcr picture we have 



^ = - l -[H 2 ,p 2 ]+T l D l {p 2 } + Y 2 D 2 { P2 } + 0{e- 
at h 



'-At\ 



(38) 



where the Hamiltonian, augmented by the Lamb shifts, 
and the two dissipation channels read: 



H 2 = —n 2 + xin 2 + x 2 n 2 {n 2 - 1), 



Di{ P } 



b 2 pb\ 



-n 2 p- ^pn 2 , 



(39) 



(40) 



D 2 {p} = b 2 p(b\) 2 \ib\fblp \ P {b\fbl. (41) 

Finally, let us gather together the conditions used in 
derivation of the reduced equation (|38| . We have 
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•x« r «~w- #«'■ < 42 > 



The validity conditions (|42j) can be recast in terms of 
the characteristic tunneling times and the nonlinear time. 
Defining T13 = h/J, t\ 2 = %/ ' J\ 2 and r„; = h/U, we 
have: ^, ^±2. < (Tt^)- 1 < 1. The rates of the two 
dissipation channels of mode b 2 have the following orders: 
Ti ~ (r 13 /r 12 ) 2 /r and T 2 ~ (t iz / T nl ) 2 /T . 

We note that a master equation similar to Eq. (|38|) 
has already appeared before in connection with one- and 
two-photon absorption in quantum optics, where its spe- 
cial cases were studied |22j. It was shown that two- 
particle absorption has properties drastically different 
from the single-particle one. In particular, the decay 



is non-exponential and, irrespectively of the number of 
particles in the initial state of the mode b 2 , number of 
particles in this mode drops to the single-particle level 
during the same time-interval [221 ] . 

Eq. (|38|) has a number of specific features. First, we 
see that in the symmetric potential (when Ti = 0) the 
decay occurs due to loss of two particles at once and 
the quantum parity, being average of the quantum parity 
operator 



P = (-l 



Mb 2 



(43) 



remains constant. For example, for the state with the 
(P(0)) = -1 (odd parity), one will have (P(t)) = -1; 
the superposition state with only the odd (even) num- 
ber of atoms will remain the state with the odd (even) 
number of atoms during all the evolution time. Second, 
for a biased potential (Fi ^ 0) there is a resonance be- 
tween the two different dissipation channels, under the 
condition Ti = 2r2, resulting in a polynomial decay of 
population. To sec this, consider evolution of the average 
population 



d(n 2 ) 
dt 



-(Ti - 2r 2 )(n 2 ) - 2r 2 (n 2 ) 2 - 2r 2 (An 2 ) 2 , (44) 



where (An 2 ) 2 = (n?,) — (n 2 ) 2 - The initial state of mode 
b 2 , to be used in Eq. (|3"5|) . is a Fock state with a good ap- 
proximation (mode 61 is emptied on a much faster time 
scale). Hence, discarding An 2 (which is justified by nu- 
merical simulations, see Fig. [3|), we get an approximation 
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2r 2 



ri-2r 2 



(n 2 (to))[l-e-( r i-2r 2 )(t-t )] 



(45) 



giving, for I?i = 2r 2 , the t _1 -decay: (n 2 (t))~ l « 
(n 2 (to)) _1 +T±(t — to). Thus, we have a quantum reso- 
nance between two different (linear and non-linear) dis- 
sipation channels of a subsystem (mode b 2 ). The match- 
ing between them is expressed in terms of matching be- 
tween the linear bias and non-linear interaction coeffi- 
cient: U = v2Ji2 = (e 2 — £i)/v2. In Fig. [3] we show 
an excellent agreement of the analytical approximation, 
Eq. (1431) . with the numerical simulations, i.e. the Monte 
Carlo wave- function method [2fJ , of the single- mode 
and the two-mode (|29|) master equations. 



V. THE MEAN-FIELD APPROXIMATION 

The quantum derivations of the reduced master equa- 
tions from the first principles, presented above, are in- 
volved. One then may inquire if the mean-field approx- 
imation, commonly applied to the many-boson models 
with large number of bosons and which is much simpler 
to analyze, can substitute the quantum derivation. Here 
we remind that the two-site model with a local dissipa- 
tion is perfectly described by the mean-field approxima- 
tion [13, [jj|. This, however, turns out to be not the case 
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A. The first reduction: equations for fa and 



Consider the strongly dissipated case of section IIIII 
The small parameter is e = fe with the conditions 
^ = 0(e) and ^ = 0(e). For t » 1/T, fa = 0(e) 
and one can integrate Eq. (|48cp by rewriting it in the 
integral form and neglecting the higher-order nonlinear 
term. Using the integration by parts, we get 
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FIG. 3. (Color online) Comparison of the two-mode mas- 
ter equation (|29[) with the single-mode master equation (|38[) 
and the analytical approximation of Eq. (|45|l . The solid 
lines give the two-mode equation, the dashed gives the single- 
mode, and the dotted - the analytical approximation. Here: 
N — 40, r_Rri2 = 200, and, from the top to bottom, 
A = UN/2J12 = 20, N/y/2, 40. To compare the results in the 
domain of their validity, i.e. after the atoms are removed from 
mode fci, we have used the Fock initial state with all atoms 
occupying mode 62 (a Gaussian with a small dispersion leads 
to the same results). 



for the three-site Bosc-Hubbard model, as we will show 
below. 

The mean-field Hamiltonian can be obtained from 



^(t)+0(e 2 ) + 0(e- — ). 

(49) 

This result corresponds to the expression for the full den- 
sity matrix (|10[) with Eqs. ([TT]) and (fT8")l . where the am- 
plitude fa is locked to that of fa and fa with the same co- 
efficient as in the full quantum case («i of Eq. (fT5|) ). Eq. 
(|49| can be now inserted into Eqs. (|48a[) and (|48b[) . We 
obtain a reduced system describing the coherent modes: 
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Hamiltonian ([3]) by replacing the boson operators with 
c- numbers (bj — > fa) Q, we get 



H = -J 13 (fa 1 fa + fafa) 

U 



Ju(fii02 + 020i) + e\Pa\ 



U\t3 3 \ 4 + -z- ([|/3i| 2 + |/3 2 | 2 ] 2 + [PI fa + /? 2 */3i] 2 ) • (46) dition on the nonlincarity 



where T^ is given by Eq. (f20f and the Lamb shift er by 
Eq. (f2"T|) (i.e., by the corresponding quantum results). 



B. The second reduction: equation for fa 

Now, let us perform the second reduction to an equa- 
tion for the amplitude fa , similar as in the quantum case 
of section IIVI In section IIVI we have assumed that the 
new small parameter is e' = 4p- with the additional con- 

^f- = 0(e')- However, let us 
for a while broaden the derivation and discard the condi- 



Note that the total number of atoms is given as 
Ylj=i \Pj\ 2 = N. The local atom loss (dissipation) part of 
Eq. (TT]) can be simply added to the mean-field Hamilto- 
nian equations as an atom loss of mode fa (as we will 
see shortly, it is describable classically; see also Refs. 
[Il[l9J]), that is 



tion on the nonlinearity. For times t ^$> 1/Tr boson mode 
fa is practically empty, i.e. fa = 0{t'). This allows us 
to simplify Eqs. (|50a[) and (|50b[) as follows 
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Thus, the mean-field equations read: 
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where, for simplicity, we have dropped the term C(e 3 ). 
Already from this system it is clear that the mean-field 
approach will not account for the boson-pair dissipation 
channel (|4ip of mode 62 , present in the full quantum Eq. 
(|38|) . Indeed, Eq. (|51a|) must be somehow integrated 
with the result to be inserted into Eq. (|51bl) . However, 
one can notice that the parameter U enters Eq. (|51ap 
only in the conjunction with a factor ~ \fa\ 2 , hence all 



the terms in Eqs. (|51a|) and (|51b|) scale as a//V (the 
scale of /3) if we fix the nonlinearity parameter A = UN, 
whereas in the quantum case the iV-scaling of the boson- 
pair channel is 0(1). 

Now, under the condition on the nonlinearity as in the 
full quantum case, i.e. S^ 2 ' = 0(e'), one can proceed 
to derive the reduced equation for the amplitude fi%. To 
this end, the system of equations for (3\ and /3* in the 
required order can be written in the following form 
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(53) 



with 7 = = ( 1 + %* 



=p— i . Eq. (f5U|) can be put into the 

integral form and integrated for times t ^$> 1/Tr, in this 
case the matrix ^M(t) enters the exponent under the 
integral, with the result 



« =-SM4 +*— > + ^ 2 >- 



We need only the first row of the matrix M 1 : 



(54) 



(M- 1 )^ = -j- 1 + 1- 2 ^\H 2 + 0{e'\ (55) 
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- 2 ^(f3 2 ) 2 +0(e' 2 ). (56) 



From Eqs. (|54 |) -([56 |) we obtain 

/3i = 7-^A +0(e' 2 ) + 0(e-^). (57) 

nT R 



Inserting this expression into Eq. (|51b[) we arrive at the 
reduced aquation for the amplitude ,02, the mean-field 
analog of the coherent n- phase mode 62: 



dfh 
dt 



£ + £)*->■* 



+0(e' 3 ) + 0(e- J ^), 



(58) 



with Ti and n\ given in Eq. (|57|. Observe that, while the 
single-boson dissipation channel, Eq. (|4T))) , is accounted 
by the mean-field Eq. (|58p (the first term on the right 
hand side), the boson-pair channel, Eq. (|4Tj) . is not. 

In conclusion of this section, we have shown that, 
while the single-boson dissipation channel of the coher- 
ent mode 62 can be described by the mean-field approach, 
the boson-pair dissipation channel cannot be captured by 
the mean-field approximation and, thus, it has quantum 
nature. 



VI. DISCUSSION OF THE RESULTS 

We have considered the derivation of the reduced mas- 
ter equations in the limit of strong dissipation on the 
example of the Bose-Hubbard model with a local exter- 
nal dissipation (i.e., the atom loss from the central site). 
The method we have used is not based on the assump- 
tion of the factorization of the full density matrix, in- 
stead we demonstrate that one can effectively solve the 
master equation directly in the lowest orders of a small 
parameter (inversely proportional to the local dissipation 
rate). On this way, one is able to obtain the reduced mas- 
ter equations for the subsystems (the coherent modes) of 
the higher-order in the small parameter (e.g., we have 
derived the equation up to the third order). 

The derivation reveals the following features. First of 
all, the full density matrix is not factorized (which is 
the usual assumption, sec for instance Ref. [131) but is 
expressed in the form of a "dressed" factorized density 
matrix, where the population of the strongly dissipated 
mode depends on that of the other modes. Neverthe- 
less, the reduced density matrix of a subsystem is shown 
to satisfy a master equation in the Lindblad form. This 
feature appears in the two reduced master equations de- 
rived in the present paper, thus suggesting an universal- 
ity. Moreover, the Lamb shifts and the dissipation rates 
of the subsystem turn out to be given by the similar 
expressions in the two cases, suggesting even the univer- 
sality of the form of expressions for these quantities. 

We have analyzed the relation between the full quan- 
tum derivation of the reduced master equation for the 
density matrix of a subsystem and the mean-field adi- 
abatic elimination procedure. We have found that the 
mean-field approximation applied to the Bose-Hubbard 
model cannot capture all dissipation channels of a subsys- 
tem, even if the external dissipation applied to the system 
is describable classically (in our case, the single-boson 
dissipation channel describing the removal of atoms). 
Namely, in the three-site Bose-Hubbard model the ir- 
phase coherent mode has the boson-pair dissipation chan- 
nel, which is not captured by the mean-field approxima- 
tion, and the single-boson dissipation, captured by it. 
This is a quite distinct situation from the two-site model, 
where the dissipation dynamics is described by the mean- 
field approximation with a good accuracy [H, Gjl • Here 
we note, however, that in the two-site model with dissi- 
pation the boson-pair tunneling is not possible. Though 
we have considered the three-site Bose-Hubbard model, 
the boson-pair dissipation channel is a general feature of 
the multi-site models, since it is the result of a nonlinear- 
ity of the model and the fact that more than one filled 
site is coupled to the dissipated site, which serves as a 
common reservoir. 

The failure of the mean-field approximation to the 
quantum master equation means that the dissipation 
channels not accounted by the mean-field have a genuine 
quantum nature. This imposes severe limitations on the 
applicability of the semiclassical mean-field approach to 



the Bosc-Hubbard models considered as open quantum 
systems, even when the external dissipation is describ- 
able classically. Invoking the link to the discrete WKB 
approach, mentioned in the Introduction, one has to de- 
velop a more general approach by going to the first-order 
approximation in the effective Planck constant 1/-/V, i.e. 
a 1/iV-correction to the mean-filed equations, to capture 



the boson-pair dissipation channel. This is an important 
conclusion, in view of the current experiments on the 
single-site addressability and controlled measurement via 
the local atom losses of Bose-Einstein condensate loaded 
in the optical lattices [12, LL3( > describable by the open 
multi-site Bosc-Hubbard models. 
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